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Abstract 



The exact solution for the shape and gravitational field of a rotating two-layer Maclaurin ellipsoid of 
revolution is compared with predictions of the theory of figures up to third order in the small rotational 
parameter of the theory of figures. An explicit formula is derived for the external gravitational coefficient 
\l J2 of the exact solution. A new approach to the evaluation of the theory of figures based on numerical 

o 

^+ integration of ordinary differential equations is presented. The classical Radau-Darwin formula is found 

f~^ not to be valid for the rotational parameter £2 = n^/(27rGp2) > 0.17 since the formula then predicts 

o 

I a surface eccentricity that is smaller than the eccentricity of the core-envelope boundary. Interface 

" eccentricity must be smaller than surface eccentricity. In the formula for £2, ^ is the angular velocity 

"^ of the two-layer body, p2 is the density of the outer layer, and G is the gravitational constant. For 

an envelope density of 3000 kg m^'^ the failure of the Radau-Darwin formula corresponds to a rotation 



X 

Cu period of about 3 hr. Application of the exact solution and the theory of figures is made to models 



of Earth, Mars, Uranus, and Neptune. The two-layer model with constant densities in the layers can 
provide realistic approximations to terrestrial planets and icy outer planet satellites. The two-layer model 
needs to be generalized to allow for a continuous envelope (outer layer) radial density profile in order to 
realistically model a gas or ice giant planet. 
* Retiree 



1 Introduction 



Kong et al. ( 2010 1 have presented an exact theory for the rotational distortion of a rotating two-layer spherical 



body with a constant density core surrounded by an envelope (outer layer) with a different constant density. 
The solution for the case when the core and envelope have equal densities, the Maclaurin ellipsoids, was 
obtained more than 250 years ago and attracted the attention of such notables as d'Alembert, Clairaut, Euler, 



Laplace, Legendre, Poisson and Gauss. The solutions are discussed in Chandrasekhar (1969 Ellipsoidal 



Figures of Equilibrium) and Lamb (1932 Hydrodynamics). Until the work of Kong et al. (2010), the classical 



solution for the constant density Maclaurin ellipsoid had not been generalized to a body with non-uniform 
density. Instead, approximate solutions, for bodies with density increasing from the surface to the center, 
have been developed by geophysicists interested in the internal structures of the Earth and planets. The 



approximate solutions fall under the umbrella of the theory of figures ( Zharkov and Trubitsyn 1978 Physics 



of Planetary Interiors) and rely on the smallness of a rotational parameter that measures the distortion 



of a rotating body. In this paper we refer to the theory developed by Zharkov and Trubitsyn (1978) for 



the shapes of fluid bodies as the theory of figures. The rotational distortions of the Earth and planets are 
indeed small, but it is important to accurately determine the small distortions to correctly infer the interior 
structures of the bodies. Accordingly, the theory of figures has been developed to high order in the small 
rotational parameter. 

While the two-layer spherical body with constant core and envelope densities is too simple a model for 
the careful study of most planets, the exact solution for the rotational distortion of such a body provides a 
benchmark against which the accuracy and range of validity of approximate solutions and numerical models 
can be evaluated. In this paper we compare results from the exact solution for the rotating two-layer 
Maclaurin ellipsoids to those obtained using the Radau-Darwin approximation and other simplifications 
based on the theory of figures. We derive formulae for the rotationally distorted two-layer sphere using the 
theory of figures valid to third order in the small rotational parameter and assess these results against the 
exact solution. Application of the exact theory is made to the Earth and planets while keeping in mind the 
limitations of using a two-layer model to represent these bodies. Both the theory of figures and the exact 



theory of Kong et al. (2010) assume hydrostatic equilibrium of the rotating body. 



2 Exact Solution for a Rotating Two-Layer Spherical Body 

We consider the distortion of a two-layer spherical body rotating with constant angular velocity fi. The core 
has radius ri and constant density pi. The core is surrounded by a spherical shell envelope of outer radius 



r2 and constant density p2- The exact solution for the distortion is presented in Kong et al. (2010). Here, 



we summarize only a few things necessary for making use of the theory. The problem is completely specified 



by three dimensionless parameters, the core-envelope density ratio pi/ P2, the fractional volume of the core 
Qv = (^1/^2) , and the rotation parameter 

where G is the gravitational constant. Among the quantities derivable from the solution are the eccentricities 
of the total (gravity) equipotential surfaces, and, in particular, the eccentricities of the core-envelope interface 
and the surface, Ei and E2 , respectively. The total or gravity potential is the sum of the gravitational 
potential and the rotational potential. Surfaces of constant total (gravity) potential are shown in Figure fl] 
for the cases (a) Qv = 0.5, pi/ P2 = 2, £2 = 0.18 and (b) Qv ~ 0.25, pi/ P2 = 2, £2 = 0.05. The eccentricity 
of these total potential isosurfaces is plotted as a function of radius in Figure [2] The eccentricity of total 
potential isosurfaces generally decreases inward except for the region of the interface where eccentricity 
changes rapidly and non-monotonically. Eccentricity has a local maximum at the interface and a local 
minimum in the envelope just above the interface. Eccentricity decreases monotonically with decreasing 
radius inside the core; the decrease is very gradual until the center of the core is approached. 

The coefficient J2 of the external gravitational field is an important quantity to determine because it 
can be measured by a spacecraft flying by or orbiting a planet. For the rotating two-layer body, J2 can be 
found from the exact solution by proceeding as follows. In a spherical coordinate system, the axisymmetric 
gravitational potential V^ outside a uniformly rotating body can be expanded as 
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where M is the mass of the body, P2n are the Legendre polynomials, R^ is the equatorial radius of the body, 
r is the radial distance from the center of the body (r > i?c), and 6 is the colatitude of the observation point 
with respect to the rotation axis. On the spherical surface r = R^, the expansion becomes 

V^{r = R,,e) - -^ [1 - J2P2(cos0) + •••]. (3) 

By determining Vg{r — Rc,0) from the exact solution, J2 can be calculated from the projection of the 
potential onto the expansion ([3]). 

The gravitational potential expansion in the spheroidal coordinate system employed in the exact solution 
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(Kong et al. 2010). In Q, ^, 77 are spheroidal coordinates with focal length C2, Pi and Qi are Legendre 
functions of the first and second kind, respectively, and i is the square root of —1. In addition p' is the 
density as a function of ^', 77', and ^o is the value of ^ at the outer free surface. Evidently, we need to make 
the transformation from spheroidal coordinates to spherical coordinates before computing Vg{r,9) from (|4|. 



The relationship between the spherical and spheroidal coordinates is (Kong et al. 2010) 

rcose = C2V(1 + C2)(1-^'), 
rsin^ — C2^?7 



(5) 
(6) 



Taking r ~ R^. and using the fact that C2 = RcE2, we find the transformation in the form 

cos 9 = E2^J{l+e){l-v''), (7) 

sine = ^2^?? (8) 

Equations (It]) and (Is]) enable us to derive 1^ and 77 as functions of 9 
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With ([9]) and (10), we are able to derive the gravitational potential as a function of the spherical coordinate 



When Vg{r — Re, 9) is available, we project it onto the spherical harmonic expansion to obtain J2. We 
expect the expansion in the form 



T/g (r = R,,9) = Y. \l'^QPi{cos9), 



where 
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By comparing (11) and (IS]), we find 



from which J2 can be simply calculated as 



1 _ GM 

5 _GM 
2^' IT " 



J2 = 



-5C2 
Co 



(11) 
(12) 

(13) 
(14) 

(15) 



3 Comparison with the Radau-Darwin Approximation 

The Radau-Darwin approximate formula can be used to predict the flattening or eccentricity of the outer 



surface of a rotating body (Radau 1885 Darwin 1900). The formula relates the normalized axial moment 
of inertia C /Ma\ {C is the axial moment of inertia around the rotation axis, M is the total mass of the 



body, and 02 is the equatorial radius) to the second degree Love number /i2 
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( Zharkov and Trubitsyn[ 1978). The Love number gives the flattening of the surface /2 = (02 — C2) /a2 by 

qh2 



/2 = 



(Zharkov and Trubitsyn 1978), where q is the small rotational parameter 
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The flattening and eccentricity of the surface of the body are related by 

h = l-{l-ElY" (19) 

Equations ([16)-19) are often used in planetary physics to determine the moment of inertia of a body 



whose shape (flattening), rotation rate, mass, and equatorial radius are known. With /2 and q known, (17) 
gives /i2 and (16) gives C/Ma1. Alternative to the flattening, the gravitational coefficient J2 can be used to 



infer the axial moment of inertia since to a first approximation J2 and /2 are related by 



(20) 



(Zharkov and Trubitsyn[ 1978), where m is the small rotational parameter given by 



GM 



(21) 



and S2 is the mean radius of the body (the radius of a spherical body with the same volume as the rotationally- 



distorted body). The Radau-Darwin formula can be rewritten in terms of J2 and m by using (17) and (21 ) 
to elinunate /i2. 

The Radau-Darwin formula can be used to predict the fiattcning or eccentricity of the surface of the 
rotating two-layer spherical body and the result compared with the value of the eccentricity from the exact 



solution of Kong et al. (2010). The normalized moment of inertia of a two-layer sphere is 
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This can be rewritten in terms of the dimensionless variables that characterize the solution of Kong et al. 
C 



(2010) as 
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Given Qv and pi/p2i the dimensionless moment of inertia is calculated from (23), and /12 follows from (16 1 



The flattening and eccentricity of the model surface is then obtained from (17) and (19 1 



Table [1] compares the exact solution for the eccentricities of the interface Ei and surface E2 (evaluated 



using the theory of Kong et al. (2010)) with the eccentricity of the surface from the Radau-Darwin formula 



E2 for the case Qv = 0.5 and pxjpi = 2 for different values of the rotation parameter £2- The agreement 
between the approximate and exact solutions is quite good in this case for all the rotation rates considered. 
Strictly speaking, the Radau-Darwin approximation can be said to be invalid for values of e2 > 0.17 since 
the Radau-Darwin formula predicts surface eccentricities less than the interface eccentricities, when in fact, 
the interface eccentricity should be less than the surface eccentricity. That this is the case can be quali- 
tatively understood by considering the flattening of a rotating sphere of uniform density p. The flattening 
is proportional to p^^. In a body that has density increasing with depth, the flattening or eccentricity of 
equipotcntial surfaces should accordingly decrease with depth. 

4 Comparison with the Theory of Figures for the GeneraUzed 
Roche Model 



The generalized Roche model is a special case of the two-layer model of this paper in which the envelope 
density p2 = . An analytic formula for the flattening of total potential isosurfaces, correct to second order 
in the rotational parameter m, is given in Zharkov and Trubitsyn ( 1978[ ). (We also derive this formula later 
in Section 5.8.1, which discusses the generalized Roche model from our perspective on the theory of figures.) 
The flattening of the surface according to this formula /J°^ is given by 
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where Pc = '"i/'"2 = Qv • The rotation parameter m is defined in (21)for S2 = r2. It is related to the 
parameters of the exact theory by 



3(^2^2) 



(25) 



2QvPi 
where 62/02 is independent oi p2 (see (IT])). The flattening of the interface according to the second order theory 

of figures ToF is 
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Equation (26) follows from (58) and (87 1 and (91) with P = (3c in the expressions for Fi and F2 



Table |2] compares the eccentricity of the surface and the interface of several models computed using the 



exact theory with values obtained from (24) and (26 1 for the generalized Roche model with the help of (19) 
to convert to the eccentricity -Ej"^ . For the exact theory we consider the sequence of values of P2I P\ equal 
to 10^^, 10^'', and 10^"* to compare with the Roche model result that takes p2 = 0. For the cases considered 
in the table the theory of figures to second order in m only slightly underestimates the surface and interface 
eccentricities. 



A formula for the surface flattening correct to third order in m is given by combining / (/3c) from (26) 



with _F3i given by (92) 



5 First-Order Theory of Figures for Synchronous Rotation and 
Tides 

The theory of figures dates back to Clairaut, who in 1743 derived an intcgrodifferential equation for the 
flattening of a rotating body in hydrostatic equilibrium (HE), but with a non-uniform density distribution 



in the interior (Kaula 1968). Clairaut's theory represents a first order perturbation theory to a non-rotating 
spherical configuration with arbitrary density distribution in layers, the level surfaces, or the surfaces of 
constant total potential, the sum of the gravitational potential and the non-inertial centrifugal potential 
(zero in the spherical configuration). The small rotational parameter m for the perturbation theory is 
related to the body's rotation period P, its mean radius i?, and its total mass M by 



ra = 



¥) (^) (--P--@) 



(27) 



where G is the gravitational constant given by (6.674215 ± 0.000092) x 10 ^^ m^ kg s ^ (Gundlach and 



Merkowitz 20001. Often, the combination GAf is known from orbital dynamics to greater accuracy than G 



itself, and the total mass is a derived parameter given by GAf /G, accurate to essentially the same fractional 



error as G, or 14 parts per million. The mean radius R (also denoted by S2 in (21 )) is the radius of a uniform 
sphere equal in density to the planet's mean density po or, 

^° - 4^aR^ ^^^> 

The flattening / is defined in terms of the equatorial radius a and polar radius c by, 

/ - — (29) 

a 

In the first order Clairaut theory, the rotating mass configuration consists of continuous layers of concentric 
ellipsoids of revolution, each with its own density and fiattening. The configuration is defined in terms of 
a single variable, the normalized mean radius (3 — s/R, which labels the level surfaces. The actual mean 

7 



radius s, in metric units, can always be recovered by multiplying j3 by the given mean radius R of the planet. 
The density p{P) of a particular level surface can also be normalized to the given mean density, so that 
5{13) = p(/3)/po. 

For perturbations of higher order than the first, these definitions can all be retained. However, the level 
surfaces are no longer ellipsoids of revolution. Even so, the radius r of a particular level surface can be 
expressed in terms of the mean radius s by a distortion to the usual polar coordinate equation for an ellipse. 
For the third order theory, the radius r can be expressed as a function of the polar angle 9 (colatitude) , the 



flattening /, and two higher-order, spheroidal shape parameters k and /i, as follows (Zharkov and Trubitsyn 

3 



1978). 
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It is convenient to express the radius in terms of /^, the cosine of 9. Then the mean radius is given by 
the integral 



/ = 



1 



r^ (^) d^ 



and the expression for the radius becomes 
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Consistent with (30) and (31), the coefficients to third order are 
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This is an alternating series that converges absolutely, and the error in the partial sum is less than the 
absolute value of the next term in the series. 



Equation (32) is fundamental to the theory of figures. With it, a coordinate transformation between 
Cartesian coordinates {x,y,z) and generalized coordinates {l3,(J3,fi) can be defined by 

X = rcos(f)\/l — y? 



r sm 



r/i 



0^1 -M' 



(34) 



Because of axial symmetry, the azimuthal coordinate (j) can be integrated out of the problem immediately. 
The series coefficients a2i are functions of j3 only, by means of the shape functions f,k,h. The differential 
volume element dr can be found from the Jacobian determinant of the transformation. 



For example, after integrating for the coordinates (j) and fj,, the gravitational coefScients J„ in the external 



gravitational potential are given by Zharkov and Trubitsyn (19781, 



Jn^- 5[P) P'^Pn M dr = - / Sil3)d [/3"+3</)„ (/3)] 



(35) 
where P„ is the Legendre polynomial of degree n. The functions 0„ are derived in IX] to third order by means 



of a definite integral in /i from minus one to plus one. The only integral that remains in (35 1 for further 
evaluation is the integral in /3, which depends on the given density distribution S{/3). 

Similarly, the principal moments of inertia, C along the polar axis, and A along an equatorial axis, can 
be expressed as third-order series. All quantities are thereby normalized to the total mass M and powers of 
the mean radius i?, so that the external gravitational potential function V is expressed as 

2i 



V = 



GM 
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(36) 



The measured gravitational coefficients J2i can be determined from orbital dynamics, as is GM. They are 
most often referred to a reference value for the equatorial radius a,.cf . However, they can be referred to R by 



multiplying each observed value of degree 2i by (orcf/i?) , consistent with the computed values from (35) 



The usefulness of the theory of figures is that a reasonable density distribution d{(5) can be found that 
minimizes the weighted sum of squares WSOS for the measured coefficients, the method of weighted least 
squares. The minimization function in terms of the observed coefficients J2i, along with their standard errors 



CT2i , and the computed values J2i from (351 is 



WSOS = J2 



i=l 
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(37) 



If the orbital dynamics is limited, such that there are strong correlations among the coefficients, the mini- 
mization can be generalized to include the covariance matrix T for the J2i from the analysis of the orbital 
data. The residuals J2i — J2i are placed in a column matrix z and WSOS is defined by the matrix operation 
z^r~^z, with z^ the transpose of z. 

In principle the theory can be extended to arbitrary order in the small rotational parameter m. However, 
it becomes quite cumbersome for orders greater than three. Nevertheless, all expressions required for a 



fifth-order theory have been published by Zharkov and Trubitsyn ( 1976 ) . In the fifth-order theory the radius 



r is still normalized to the mean radius s, but (32) is replaced by the following Legendre series, valid to 
arbitrary order n 

n n 

r . — . . — - 

(38) 



1 + ^ sojin^ + ^ ■m^S2jP2] (m) 
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This expression can be substituted into ( 31 ) and expanded in a power series in m to order n. Then coefficients 



of rn* for powers of i greater than or equal to 2 can be set to zero. This yields n — 1 equations in the n — 1 



coefRcients sqj, which can all be evaluated. For example, at fifth order the values are 

1 2 
S02 = -^S2 

2 3 
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504 = -TTTT (35S4 + I8S2S4) 

505 = -Y^ (33^2 + 250s2s2) (39) 



A substitution of the sqj values so determined into (38) yields the n order expression for r/s, comparable 



to (32 1, the third-order expression in the spheroidal functions /, k, and h. The parameter S2 in (38) and 



(39) is not the mean radius parameter used in (21 ) 



The easiest way to compare the theory with observation is to refer the measured gravitational harmonics 



to the mean radius of the planet. The advantage is that (35) directly represents what is being measured. 
However, the mean radius of a planet depends on its rotation period, which is not always known, and which 
may not even be constant throughout the interior. In this sense, the measured equatorial radius is a more 
fundamental observational constraint. Therefore, when the measured J„ are referred to the equatorial radius. 



the theoretical values given by (35) must be multiplied by the ratio (s/a)'^. This ratio can be found from 
the inverse of the expression for r/s, with fi set equal to zero. With this approach, care is required in order 
to make sure that higher order terms in (s/a)" do not enter into the theoretical expression for J„ and bias 
it. Furthermore, this is true in general. If one is computing a first-order Clairaut spheroid, it is important 
to make sure the series are always strictly truncated at first order. The same can be said for a second-order 
Darwin spheroid, or to any spheroid of arbitrary order n. 

5.1 The Level Surfaces 



The application of (35) requires that the functions S2j{l3), or alternatively /(/3), h{f3), and A;(/3) in the 
expression for r/s of ( |32[ ), be found for a given interior density distribution, expressed in terms of the 
normalized density (5(/3). This is accomplished by finding the level surfaces on which the interior gravitational 
potential is a constant. The interior potential at a normalized mean radius /3 is determined by the amount 
of mass interior to the level surface, an integral over the volume from zero to /3, plus the amount of mass 
exterior to the level surface, an integral from /3 to one. Let the gravitational moments of the mass lying 
internal to /3 be given by S2j{/3). The normalized potential Vq from the mass interior to /3 is given by the 
expansion of s/r in Legendre polynomials. We illustrate this procedure for the spheroidal functions to order 
three. The more general procedure for the S2j functions to arbitrary order is similar. 



The first step is to invert the expression for r/s in (32) and expand it in a series to order three in the 



small rotational parameter m. The result is a function in even powers of fi. Next, the Legendre polynomials 
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in /i can be inverted to any arbitrary degree to obtain powers of /i in terms of the polynomials. For third 
order in m the result is 

m' = ^ [1 + 2P2 (p)] 

// = ^[7 + 20P2(m) + 8P4(m)] 



The next step in the procedure is to substitute the powers of fi given by (40) into the series for s/r. The 
result is an expansion of s/r in a series of Legendre polynomials in the form 



^0=-=E^2,^2,(m) (41) 

J=0 



For the spheroidal functions, the coefficients Cj,. can be written as 

^° ^ ^+45^ +2835^ +315^^ 
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^«°- 231 ^^+231'^ -231^^ (42) 

This completes the expansion for the zero-degree gravitational moment So, which is basically a mass function 
given by 



3 



/3 



^o = ik / ^''5(z)dz (43) 

P Jo 

For any interior density distribution given by 5{/3), the function Sq must be equal to one at the surface of 
the planet, where /3 is equal to one. Otherwise the interior model will not be consistent with the observed 
mass and mean radius. 

In general, the gravitational moments 6*2; are included in the level-surface theory by series expansion in 



powers of the inverted r/s in (41 1 times the appropriate higher-degree Legendre polynomial, as follows 

V,^(^'-y^^'p2^i^i) z = 0,l,2,..-,n (44) 

The series expansion to order m for a particular degree 2i is carried out to order n — i. The powers of /i 



given by (40) are substituted into the series for Vi. The result is an expansion in Legendre polynomials that 
can be written 

n 

^^ = E ^2j^2, (m) * = 0, 1, 2, • • ., n (45) 



3=0 



The evaluation of the coefficients Cj, for the spheroidal functions is given infAlfor orders 1, 2 and 3. The 



coefficients for order zero are given by ( 42 ) 
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The gravitational moments 5*2, for the potential exterior to the level surface labeled by /3 require potential 
functions VI, which are defined by 

2i 



'=(^)V2,(m) * = 0,l,2,...,n (46) 






After similar expansion in powers of m as for Vi, the potentials for mass between /3 and the surface at /3 = 1 
can be expressed in terms of coefficients Cj' by 

n 

^' = E ^2',^2, {^i) i = 0, 1, 2, • • •, n (47) 

The coefficients d^j are given to third order infAJ 

So far we have been concerned with the expansion of the internal gravitational potential to order n in a 
series of Legendre polynomials of degree 2n in the general coordinate ^. However, the planet is in rotation 
about its principal axis of maximum moment of inertia, the z axis. In this rotating non-incrtial coordinate 
system the planet deviates from a sphere because of a centrifugal force generated by a rotation in inertial 
space, a rotation with respect to the "fixed stars" . Relativistic corrections to this Newtonian model are 
ignored in the theory of figures for planets. Therefore the centrifugal force per unit mass can be represented 
by the following potential function V^ot 

V;ot = ^^)%^sin2 (48) 

This potential function can be made consistent with the gravitational potentials Vi and VI by replacing the 



period P by the smallness parameter m according to (27), by replacing ^ by the Legendre polynomial P2 



according to (40), and by normalizing to the gravitational potential GM/R at the surface. The result is 



(ZharkovandTrubitsyn 1978) 



Q-^m(^)'[l-P2(/i)] (49) 



The centrifugal potential Q can be expanded to arbitrary order by means of ( 38 ) or to third order by (|32 



The third order coefficients corresponding to the third order coefficients for Vi and VI are 

^ \i 45'' 189-' 315 

/I 20^ 38 ^, 16, 
^^ yi 63-' 189-^ 45 

_ / 8 „ 76 „2 32 , 

^ V35 231 -^ 55 

The small rotational parameter m enters explicitly in the theory of figures by means of the centrifugal 
potential Q. 
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All the coefficients derived so far can be collected into complete expressions for the internal potentials 
Ai on the level surface labeled by /?. These potentials to arbitrary order can be written as 

n 

A^^T.i^hSy+^2JS!,,)+Q, z = 0,2,4,...,2n (51) 

3=0 



The gravitational moments can be written in terms of the following integrals (Zharkov and Trubitsyn 1978) 

S[^r' f5{z)d[z'-^c^[] (52) 

■7/3 

The functions (pi and (j)[ represent the integral of the gravitational moments over /i as follows 

<^. = ^£p.(M)Q"'dM z = 0,2,4,.,2. 

</>: = ^^^/^^.(M)(^)'"'dA* z = 0,2,4,...,2n (53) 

When i is equal to 2, the integration for (p^ must be carried out as a special limiting case. The integral is 

3 



02= 2^ ^2(/i)ln(^) d/i (54) 

The functions under the integrals for (pi and (p[ can be expanded to arbitrary order in m and integrated. 
The results to order three are given in ffl Results to order 5 by means of ( 38 ) are given by Zharkov and 



Trubitsyn (1976) 



Evaluations of Aq^ A2, A^ and Aq are given in Appendix B. In order that the potential be a constant 
on level surfaces, all potentials of order greater than zero must be zero. This means that any Ai with i 
equal to 2 or greater can be multiplied through by a constant. It also means that A2 can be used to solve 
for m, as it appears explicitly in Q2. By substituting this value of m into the higher-degree potentials, A4, 
Aq and higher, they can be simplified. They do not contain m explicitly. Aq is the only potential function 
that is not zero. For this reason it represents the total internal potential at normalized mean radius /3, with 
the centrifugal term included in the potential. It is the potential that enters in the equation of HE. The 
pressure p{l3) on a level surface and the total gravitational potential U{/3) can be normalized by the following 
relations involving the given mass M and mean radius R for the planet 

RpiP) 



xiP) 



GMpo 



In terms of these normalized variables, the equation of HE in the interior is given by [Zharkov and Trubitsyn] 



(1978) 



g-^ 
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5.2 Solution to the Level-Surface Problem 

The objective of a solution to the level-surface problem is to find the gravitational moments and the shape 
of the planet at its surface, and to compare the theoretical result with what is observed for the surface shape 
and external gravitational field. This result depends on the interior normalized density distribution S{l3), 
which can be a given function, as in the two-zone model treated here, or it can be obtained from a known 



equation of state (EOS) by including ( 56 ) in the solution of the overall problem. When the EOS is given for 
the internal material as a function x{^)j most likely in zones, the following differential equation for 6 can be 
included in the solution for the theoretical interior model 

dx\ cM _ d(/3^Ao) 

dSjdp-^ d/3 ^^^> 

This suggests that it might be useful to have the level-surface problem not in the form of integrodifferential 
equations, but in the form of differential equations only. An advantage of this approach is that a numerical 
solution to a set of differential equations (ODE) can be carried out to high precision, in fact to far more pre- 
cision than needed to justify the accuracy of the observational constraints on a static model. An alternative 



method used in a previous paper (Anderson and Schubert 2007) can cause precision problems. The method 



expresses the shape coefficients /, k and h, or S2i, as polynomials in /?, and forces the polynomial coefficients 
to satisfy the equations A2i = 0. The problem with this approach is that a finite number of polynomial 
coefficients can never be found that satisfy the equations everywhere on the interval < /3 < 1. Numerical 
compromises must be made in order to satisfy the equations on average over the interval. With the ODE 
approach, the solution for the shape parameters can be automated. 

Using this ODE approach, we first express the shape coefficients as power series in the small rotational 
parameter m. We illustrate the method for /, k and h, and use it for the two-zone model, but it can be 
extended to higher orders as well. The three spheroidal functions can be written as 

/(/3) = mi^i(/3) + m2^2(/3) + m^FsiP) 

ki(3) = m2X2(^)+"i3i^3(/3) 

ft(/3) = m^H^iP) (58) 

The first step in the procedure is to substitute these expressions for /, k and h into the functions (f)2i and 
4>2i given in Appendix C and to drop terms of order higher than three. The next step is to substitute the 



resulting power series into the expressions for 5*2^ and S'jj given by (52). Finally substitute the resulting 
gravitational moments and the shape functions /, k and h into the expressions for A2, A4 and Aq given in 
Appendix B. Then expand to order three in m. This completes the setup of the level-surface problem for 
the ODE approach. 

The lowest-order level-surface potential is A2 to the first order in m. Call it A21. It is obtained as the 
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coefficient of m for A2 from the setup of the problem. It can be written as 

^21 - -;^ + 5oFi {(3)-l [ 6 (z) F[ (z) dz - -|. / ^"'5 (^) (5^1 ('^) + ^^1 (^)) d^ = (59) 

2 5 Jfl bP'^ Jo 



where S'q can be evaluated by the integral of (43 1. Differentiate A21 once with respect to /3 to obtain 

3 f' ....._ .. . 1 



^2i = ^y z45(z)[5Fi(z)+zF{(z)]dz-^5o[3^i(/3)-/3F{(/3)]=0 (60) 

Multiply this derivative by /3^ and differentiate once again. The result is 

(/3«A^0' = 6/3'^ (/3) [Fi (/3) + /3F{ m - /3^S, [6F, {(3) - fi'F[' (/3)] = (61) 



This last equation (61) is Clairaut's differential equation for the flattening function. However, the first two 



equations contain integrals that are not known, and they will enter into higher-order ODE. Therefore, we 



solve for the two integrals from the two equations ( 59 ) and ( 60 ) and at the same time solve for the second 



derivative of Fi from the third equation (61). This establishes a procedure for all higher orders. The result 

is 

J^^ z^ 6 (z) [5i^i (z) + zF[ (z)] dz = IS^SoFi (/?) - ^/3^5oF( (/?) 
f^ S (z) F[ (z) dz - -^ + ^^oi^i (/?) + ll3SoF[ (/3) 



A.,.,-A(^),,,_|(^. 



Fi' iP)-^F,{P)--(^]F,{/3)^-(^]F[{(3) (62) 



The ODE in (62) can be solved for Fi and F{ and the result can be substituted into the two integrals. The 
solution to the ODE to first order in m and the corresponding two integrals are now available for higher 
order ODE. The boundary conditions on the solution are discussed in section |5.3| and they are applied to 
the two-zone model in section p3.6| Note that the density function that completely determines -Fi is given by 
the ratio S/Sq. 

The next function for consideration is K2. It is derived from the coefficient A42 of m'^ in A4. This 
time A42 is divided by /3^ and differentiated. Then the result of that first differentiation is multiplied by 
P^'^ and differentiated once more. This establishes the procedure for all the shape functions. When Aq is 
involved, it is divided by /3^ and differentiated. Then the result of that differentiation is multiplied by /3^^ 
and differentiated once more. The procedure can in principle be carried to higher orders. For each shape 
function, three equations are solved for two unknown integrals and the second derivative of that particular 
shape function. The sequence of steps for deriving the ODE is Fi, K2, F2, H3, K^, F3, and so forth. The 
result can be expressed as a nonlinear homogeneous differential equation plus a function of (3 that is built 
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up by means of the sequence of derivations. We express the ODE in the form 

/3'^r + 6/3('|^')F/-6('l-|-')F, = G2. z = l,2,3 
P^K': + 6/3 (^\ K[ - 2 Ao - i^\ K, = G4. z = 2, 3 

P^H'; + (Sp(^\H[-&[7-^\H, = Ge^ * = 3 (63) 



These equations differ in form because of the way k and h are defined in ( 30 ) . The functions Gji are given 
in Appendix D. 

5.3 General Boundary Conditions 

The derivatives of the shape functions at the surface where /3 is equal to one can be found sequentiahy, 
similar to the technique for finding the ODE. For the function Fi, the potential A21 is multiplied by f3^ and 
differentiated. This is done for A22 and ^423 as well. For A42 and A43 the multiplier before differentiation 
is /3^, and for Ags it is 0^ ■ The resulting derivatives are evaluated for /3 equal to one. Consequently, the 
integral with limits of integration from /3 to one is set to zero. The integral representing Sq is evaluated at 
the surface such that Sq is equal to one. The second derivatives are eliminated by substitution of the ODE, 
again evaluated at the surface. The resulting first derivatives of the potential functions multiplied by the 
appropriate /3* can be set to zero and all the derivatives of the shape functions at the surface boundary can 
be found sequentially. As a result, the surface boundary conditions are given by 



^21 = 


1~\f..~^K2. 




^21 = 


-S + S^-n^-^^-4^- 




^31 = 


f + ^^11 - 5Fi\ - 6i/3i - 2i^2i 




^31 - 


25 25 ^ 137 ^2 5 ^ 12 904 
24 + 168^" + 168^" 4^^^ + 11^^^ + 231^'^ 


524^ ,, 
-231^"^-- 


Fk- 


155 143 47 2 ,73 ,19 2 

72 84^" 147^" + 9^" + 42^^^ + 3^"^^^^ 

11^^^ +2695^"^'^ +7^^^ 


- 2F31 ~ ^^31 



AK; 



31 



(64) 

This gives the derivatives of the shape functions at the surface. One more set of boundary conditions 
is needed for a unique solution to the ODE, and hence for a unique interior model for a given density 



distribution (5(/3), or for a unique EOS that can be integrated by (57) to yield a unique density distribution 



One approach is to iterate on the surface functions, which must satisfy the boundary conditions given by 



(64), until finite functions are obtained at the origin. However, this iterative process can be tedious. An 
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alternative, which we adopt here, introduces a constant-density core into the interior density distribution. 
The core radius Pc can take on any value in the interval < /3c < 1, and in principle it can be arbitrarily 
small. However, as the core radius approaches zero, the numerical precision required to evaluate the shape 
functions at the core boundary becomes arbitrarily large. For every model, there is a practical lower limit 



to the core radius /3c. We illustrate this method of a core plus overlying envelope in Sec. 5.6 



5.4 Calculation of the Normalized Pressure in the Interior 

To the first order in ni, the pressure depends only on the density distribution. The differential equation for 



X is simply ( Zharkov and Trubitsyn 1978 ) 



dx 

d/? 






(i5{p) (65) 



The boundary condition for a solution to (65 1 is that x is equal to zero at /3 equal to one. The density 
distribution can be piecewise continuous, as in the two-zone model considered here. However, the pressure 
and gravitational potential must be continuous over a density discontinuity. This implies that the spheroidal 
functions /, k and h (or the S2i functions) and their first derivatives must be continuous throughout the 
interior. In addition, the derivative d(5/d/3 must be less than or equal to zero throughout the interior, so 
that the density either remains constant (incompressible material) or increases with depth. Also (5(/3) must 



satisfy the boundary condition that the gravitational moment Sq as given by ( 43 1 must be equal to one at the 
planet's surface. The surface is defined such that all the planetary mass is contained within the outermost 
level surface with /3 equal to one. Even so, the pressure and the density can be made to match a model 
atmosphere. The atmosphere is a part of the total mass. In that sense, it is more realistic to define the 
surface at the 100 mbar level in the atmosphere, near the top of the troposphere, not at a more standard 
one-bar level. Nevertheless, the one-bar level is an acceptable approximation to the surface. At least this 
approximation avoids the complication of treating the atmosphere as a separate zone in the level-surface 
computation. There is something to be said for separating the atmospheric modeling from the interior 
modeling, and simply making sure the two are consistent at the one-bar level. For one thing, the atmosphere 
is not static, but is dominated by observed zonal flows for all four giant planets in the solar system. The 
theory of figures is a static equilibrium theory. A level surface of one bar in the atmosphere is stretching 
the static assumption as it is. It is a reasonable level to stop the interior modeling. In order that both the 
density and the pressure go to zero at the surface, the density must go to zero at the surface. This introduces 
another constraint on the interior density distribution. A separate constraint is that the derivative of the 
density distribution at the surface is equal to the derivative in the atmosphere at the one bar level. With 
the inclusion of the constraint on So previously mentioned, this results in a total of three constraints on 
the interior density distribution. Physically, these three constraints mean that the total mass of the model 
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is equal to the measured total mass of the planet, and that the interior density distribution matches the 
density distribution in the atmospheric model at the one-bar level. 

By means of the derivation of the ODE for the shape functions in the interior, it is straightforward 
to derive the second and third order terms in the differential equation for the pressure. All the integrals 
necessary for an evaluation of Aq are available from the derivation of the ODE. The first order term in 



(65) contains only zero-order shape functions. Similarly, the second order terms in the derivative of (3 Aq 



contains only first order terms in the shape functions. The right side of (56) can be expanded in powers of 
m, and each order can be integrated separately for purposes of obtaining the normalized pressure x W) = 
Xo il3) + mxi (/?) + '71^X2 (/3) + Tn^Xa (/3) to third order. The four functions for the integrations are given by 

d {13^ Ao 



loo; 



d/3 
d(/32Aoi)_2 






d/3 3 

d/3 45^ ^^^' + ^^^^ + 45^^° ^^^' + ^^^'^1 + ^1^ 

d {P'Ao3) 



d/3 135 



13 [5F^ - 6F2 -I- 2/3FiF[ + /3 (/3Ff - 3F^)] 



~pSo{385F^ + 23ipF^F[) 

24 
-^^l3^SoF[ [2IF2 + 12i^2 + /3 (2/3Ff + 21F^ + 12X0] 

24 
^^^l3SoFi [105F2 + 6OK2 + /3 (25/3Ff + 21F^ + 12K!,)] (66) 



The pressure can be found by multiplying the four derivatives in ( 66 ) by the normalized density S (/3) and 
integrating, with the boundary condition x (1) equal to zero. 

The method described here can be applied to the simple case of a planet made up of incompressible 
material in HE. The normalized density is a constant equal to one, and the zero degree gravitational moment 
So is also a constant equal to one. The ODE simplify considerably, but that fact can be ignored, and our 
general numerical procedure can be applied to the constant-density case. As a result, the numerical solution 
to the ODE yields the result 

^ 5 / 15 925 2\ 

•' 4 V 56 1568 J 

h = (67) 

The normalized axial moment of inertia C /Ado? for this configuration is equal to 2/5, independent of m. 



A numerical integration of the four pressure functions in (66), with 5 {j3) equal to one, yields the following 
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result for the normalized pressure 

/ n9\ A 1 41 . 1235 o\ 

^-(^-^) (2-3™ -72^" -2268"^) (^^) 

In the above, the real numbers returned by the numerical integration have been replaced by nearby rational 



numbers with small denominator. This has only been done in (68 1. 

The next simplest density distribution is the linear distribution. Because it implies compressible material, 
it is a far better approximation to a real planet than the constant-density model. The normalized density S (/?) 
is equal to 4 (1 — /?). For purposes of applying our numerical procedure, we introduce a core of normalized 
radius /3c equal to 0.05. The normalized constant density in the core is equal to 3.85. The gravitational 
moment 5*0 in the envelope is equal to 4 — 3/3. In this model, there is a negligible fractional core mass equal 
to 77/160000. The numerical integration is carried out in the envelope over the interval 0.05 < f3 < 1. 

5.5 Calculation of the CoefRcients J„ in the Exterior Gravitational Potential 

The solution to the differential equations to third order in the small rotational parameter m yields the six 
shape functions Fi, F2, F3, K2, K^, H3. If good observations of the shape of the planet are available, such as 
for the Earth, these shape functions can be used directly to constrain the envelope density 6e- However, for 
the outer planets, the measured zonal gravitational coefficients J2, Ji and Jg provide far better constraints 



on 5e ■ The calculated values of the three coefficients are given by ( 35 ) . 

The procedure for finding values of the gravitational coefficients in terms of the shape functions is to first 
express the coefficients as a truncated power series in m according to 

J2^mJ2i + m^ J22 + TO^ J23 

J4 = m^ J42 + m'^Jis 

J6=mV63 (69) 

Next we recognize that the coefficients, when referenced to the equatorial radius a, are proportional to the 
gravitational moments Sn by 

5„=(^)"j„ (70) 



and where, from ( 32 ) with /i set equal to zero 



f--^/-^/'-^/'-S-i-> <") 



Substitute (70) into the expressions for the potential functions ^2,^4,^6 given respectively by (107), (109), 



(110), and evaluate at the surface. Use the truncated series of (69) for the Jn and the similar series for the 



shape functions given in (58). The functions S!^ are all zero at the surface. Expand to order three in m. 



Since all the coefficients of powers of m are zero, this process yields six equations which can be solved for 
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the six J„ functions in terms of the six shape functions from the differential equations, again evaluated at 
the surface. The result is 



^21 = 

J22- 
Jiz- 

Ji2- 
JiS.- 



1 



;Fn 



3 3 

2 1 2 2 8 

Fii — 77-Fii + 77-F2I + 7:^^21 



21 3 

11 2 
"147^"^ 

^11 - -ri^ii 



21^ 

2 2 2 

— _r9i — — r ^^ ^o^ -f- — _rsi 

21 3 3 



^il31 - ^i^21 

21 -^ 105 



7 



32 
35 



i^21 



22 2 43 4 8 



192 ^ 208 ^^ 

385 ^^ 385 ^^ 



20 o 8 o 80 160 128 

21^" + 7^" + 231^- - 231^- + Tf ^"^- 



40 
147 



i^iiif2i 



21 



K^i 



32 



2695^"^'^ - §if3i 



(72) 



5.6 The Two-Layer Model 

The two-layer model consists of a constant density core of normalized density 5c, plus an envelope of 
normalized density 5e- The envelope density can be a function of /3, or even piecewise continuous in two 
or more zones overlying the constant-density core. The two densities are connected by means of the mass 



constraint implied by ( 43 ) , and they must satisfy the following equation 

nl 



5c Pi + 3 



P'5e (/3) d/3 = 1 



(73) 



Pc 



A particular interior model is defined by the envelope density 5e and the core radius jSc- The core density 
5c is a derived constant obtained from (73 1. As the core radius approaches zero, the core density approaches 



positive infinity. However, the core mass given by 5cl3c is finite at the origin. For /3c arbitrarily small, the 
core mass can represent a point-mass core with mass greater than or equal to zero. Whatever the values for 



5E{f3) and (3c, the ODE of (63) can be solved exactly in the core, and the second set of boundary conditions 
for the envelope integration can be found at the core-envelope boundary. 



5.7 Solution to the Theory of Figures in a Constant-Density Core 



The functions needed for the ODE of (631 are Sq and 5 /So- For a constant-density core, 5*0 is simply 5c and 



the ratio 5 /So is one. It follows from (63 1 and (112) that the first-order flattening function Fi is a constant. 
It has the value Fib everywhere in the core and its derivative is zero within the core. This simplifies the 



other ODE of ( 63 ) considerably. The equation for K2 is 

P^K'^ + Qfm'^ - IAK2 = 



(74) 
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The boundary conditions on ( 74 1 are that K2 is finite at the origin and that it is equal to K2B on the 



core-envelope boundary. The solution to ( 74 1 and the boundary condition at j3 equal to j3c are 

2 



K2^K2B 



K'2B='^ 



m 



(75) 



Both the shape functions and their first derivatives are continuous at the core-envelope boundary. Therefore, 
the above boundary condition applies to both the core and the envelope at /3 equal to j3c- Similarly, the 



equation for F2 from (63) and (112) is 



fi'^F!^ + QI3FI2 = -8K2 



(76) 



with the solution 



F2—F2B + ^K2B 



1- 



/3c 



(K2b\ 



The equation for H^, is 



with solution 



^^^^ l\Pc) 



/327J3 + 6/3iJ3 - 36iJ3 = -i 






F1BK2B 



H'i — Hy,B 






4FibK2b 



Pc 



1- 



/3c 



"'--A^X^) 



'3B' 



The equation for K^ is 



13"^ K'^ + QliK'^ - 14/^3 = 



i) ^-''■--i<|:,l «- 



ith solution 



^ ^ ( liV 24 /^\ 
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Finally, the equation for Fj, is 
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(77) 



(78) 



(79) 



(80) 



(81) 



(82) 
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with solution 



F^—F-iB — FibK' 



Fopi — — 



3B 



"3S 

3288 
2695 



1BJ^2B 



1 



FibK2B 



1 

7 ^ 
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92 (IhB\ 8 /i^\ 



(83) 



5.8 Solution to the Theory of Figures in a Constant-Density Envelope 



When the normalized envelope density 5e is a constant, ( 73 1 yields the following expression for the normalized 
density in the core 

5c = 5e+ ^-^^ (84) 



1-5e 
PI 
The gravitational moment ^o in the envelope is similarly 

1-5e 



So = Se 



/33 



(85) 



A value for Se, and also Sq from (85), can be substituted into the differential equations given by (63) and 



(112), and the solution for the spheroidal functions in the envelope can be found subject to the boundary 



conditions at the surface, given in Sec. |5.3[ and the boundary conditions at normalized radius (3c given by 
expressions at the core boundary in Sec. |5.7| The normalized core radius /Sq enters through the boundary 
conditions of Sec. |5.7| only. For any finite value of Se the solutions to the ODE are complicated and lengthy, 
although solutions do exist for the two-zone model considered here. Nevertheless, beyond this idealized two- 
zone model, numerical integration is required when the envelope density varies with /3. Perhaps the exact 
solutions to the ODE for constant density are useful for purposes of checking the precision of the numerical 
integration, but they have no particular advantage to the problem of interior modeling. In practice, numerical 
integration is a useful general approach for any envelope density, including constant density. The precision 
of the numerical integration can be adjusted such that it is competitive with the exact solutions to the 
constant-density envelope, especially given the limited accuracy required by the observational constraints. 

However, the exact solution is tractable for the case where the density in the envelope is zero. This is 
the generalized Roche model considered by Zharkov and Trubitsyn ( 1978 ) and discussed in Section l4J The 
envelope contains no gravitational mass contribution to the HE, but there is an inertial contribution from the 
centrifugal potential, and hence a contribution to the surface shape of the planet. For purposes of illustrating 



the ODE approach, we solve this Roche case in section 5.8.1 and finally consider the general case of finite 
density in section [5.8. 2 1 



22 



5.8.1 The Generalized Roche Model 

For Se equal to zero, the density of material in the core is given by po/Z^c- ^^ turn, the differential equation 
for Fi in the envelope becomes 

f3^F[' - 6Fi = (86) 

The boundary conditions on the solution for Fi is that F[ — 5/2 — 2Fi at f3 equal to one, and F[ equal to 
zero at /3 equal to /3c- With these boundary conditions, -Fi in the envelope is given by 

2/35 + 3/35^ 



F = - 



4/32 



^,^3(^5.^5^) 



2/33 
1 
4 



Fu = ^(2 + 3/3S) 



FiB^l^'c (87) 

Here Fu is the value of Fi at the surface, and Fib is the value of Fi at the core boundary. All the quantities 



in (87 1 are needed for solutions to the ODE for the higher order shape functions. The value of Fi in the core 
is determined by the envelope density distribution, which in this Roche case is zero, but it is true in general 
for any envelope density distribution. For the Roche model, the density distribution in the core is simply a 



constant, given by Fib according to the core solution of Sec. 5.7 



The differential equation for the function K2 is obtained from (63) and (112), and after the envelope 



density is set to zero and the solution for Fi is inserted into 6*42, the equation for the Roche model is 
Again, with the boundary conditions from Sections [5.3| and [SJl the solution for K2 is 



^^-32 /34 


-Wh') 


^■' 16 /35 

K2B^0 





(89) 



Similarly, the differential equation for F2 is 



I3'F^' - 6F2 = -^ (4/35 + 11/3S) (90) 
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and the solution is 



F2 = 



i3h 
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(91) 



This process can be extended to third order, although the expressions for H^, K^ and F3 as a function of /3 
in the envelope become more lengthy. We list here only their values at the surface, which are 



^31 = ^(1-/9^)' (38+67/35,) 

i^3i=-^(42+49/3^ + 200;9S"424/3i,o - 200/3^,3+333/31,^ 
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The values for the gravitational coefficients can be obtained from ( 72 ) . The results are 
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(93) 



These results for the gravitational coefficients in the generalized Roche model agree with |Zharkov and| 



Trubitsyn (1978), except for J23- Total agreement is a good check on our ODE method, as the derivation 



m 
in 



Zharkov and Trubitsyn ( 1978 ) is quite different from ours. We suggest that the third order term for J2 



Zharkov and Trubitsyn (19781 contains typographical errors. For example, by setting the core radius to 



one in ( 93 ) , the case of a constant-density planet is recovered. The result is 
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(94) 
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This is correct (Zharkov and Trubitsyn 1978), and it is a good check on the ODE method. However, if Pc 
is set equal to one in (34.6) for J2 in Zharkov and Trubitsyn (1978), the resuh is J2 ~ {l/2)m— (5/28)m^ — 



(1889/4410)771'^. This is not correct. We conclude that there is agreement between our ODE method and 



the method of Zharkov and Trubitsyn, but only if the third order J2 term in Zharkov and Trubitsyn ( 1978 ) 



is brought into agreement with J23 in ( 93 ) 



Although the generalized Roche model is an idealization of a real giant planet, it illustrates the method. 
Starting with a density distribution in the envelope given by 5e{P) and a core radius /3c, the zonal grav- 
itational coefficients in the external gravitational potential can be calculated. In general, the results are 
obtained by numerical integration of the ODE, but the numerical values analogous to the six functions of 



( 93 ) can be calculated to any arbitrary precision. A comparison of the calculated values with the measured 



values is achieved by calculating the value of the small rotational parameter 777 for the planet in question, 



and then by applying ( 69 ) 



5.8.2 Model for a Finite-Density Envelope 

Even for this simple case of a finite-density envelope, the solution to the ODE can be obtained by numerical 
integration. For purposes of illustrating the method, we pick a normalized envelope density of 1/2 and a 
normalized core radius of 1/2. By ([73| the density 5c in the core is equal to 9/2, and the percentage of 



the total mass in the core [ScPq) is 9/16. This particular choice of 5e and Pc results in a fairly simple 



differential equation for the first-order function Fi . By ( 63 ) we have 



(3^ (1 + /3^) F{' + Qli^F[ - 6i^i = 



(95) 



The integration can be done numerically subject to the boundary conditions of sections |5.3| and |5.7[ which 
for Fi(/3) are 



F{(l) = --2Fi(l) 
i^{(^c)=0 



(96) 



The limits of integration are from j3c to one. After the numerical integration is complete, a value of Fi 
anywhere on the interval Pc < /3 < 1 can be found by numerical interpolation. This solution in the envelope 
can be matched to the solution in the core given in Sec. |5.7[ A plot of this particular case throughout the 
interior is shown in Fig. [3] 

Because the differential equation for K2 involves both Fi and its first derivative, it must be evaluated 



numerically by interpolating in the numerical solution to (95). Furthermore, the boundary condition at the 



surface is not known until Fi at the surface is known. Therefore, we do not write down the differential 



equation that must be integrated, but instead numerically evaluate it according to (63) on the interval 
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/?c ^ /3 < 1- The boundary conditions for this special case, with Se and /3c both equal to 1/2, are obtained 



from the expressions given in Sec. 5.3 and Sec. 5.7 and include the solution for Fi at the surface. In general, 
the boundary conditions are 



X^(/3c)=2 



K2 Wc] 



c 



^2(l)-Y^-4^i(l)-4if2(l) (97) 

and by numerical interpolation in the previous solution, Fi (1) is equal to 99060576/131853043, accurate to 
16 places past the decimal. With these boundary conditions, numerical integration yields the solution for 
K2 in the envelope, which can be matched to the core solution and plotted. The result is shown in Fig. |4] 



The procedure is similar for the function F2, and the differential equation from (63) involves the previous 



solution for both Fi and K2 and their first derivatives. The boundary conditions are 

F^ (1) = - A + l^Fi (1) + Ifi (1)^ + ^K2 (1) - 2F2 (1) (98) 

where for this special case, K2 (1/2) is equal to 143636/109713139 and K2 (1) is equal to 6168175/61992373, 
again accurate to 16 places past the decimal. After numerical integration, the solution for F2 is represented 
byFig.[5| 

The above process can be repeated for the third-order functions H3, K3 and F3, in that order. As each 
function is introduced, all previous solutions are used in both the ODE and in the boundary conditions. 
The results for the special case considered here are represented by Fig. |6j Fig. [7] and Fig. [8j All six plotted 
functions can be evaluated at the surface. As a result, the shape of the surface is given by ( [32| with (3 equal 
to one and with 

1424483 8128 2 14522 3 

•'^"1896036'"^ 132573"^ "^ 96145"^ 

9750 2 139 c, 
^^=97991™ -46715^" 

These expressions for /, fc, and h at the surface are accurate to 10 significant digits. For all practical purposes 
they are limited only by the uncertainty in the small rotational parameter m, and of course by truncation of 



the series at order m . Similarly, by (72 1, the surface conditions can be used to calculate the gravitational 



coefRcients in the external potential. The results for this special case are 



156041 980 , 339 

'2- 



J2— „„^ rn — ———-rn? + -———-m^ 



931420 25913 46459 
6381 n 2864 



2 I ^""^ ™3 



'^^ 56362"^ ■ 63535" 
^^27^ 3 

^ 46804 ^ ' 
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6 Application to Planets 

In this section we apply the two-layer model of this paper (with constant core and envelope densities) to 
the planets. The model is a simple one for planets, and it is most applicable to terrestrial planets and icy 
satellites that have a differentiated structure consisting of either a metallic core and a rocky mantle or a 
rock, metal core and an icy mantle. Application to the gas and ice giant planets will also be made, though a 
constant-density envelope is not a very realistic model of these bodies. However, with some generalization of 
the two-layer model to include envelopes with arbitrary radial density profiles, application to giant planets 
can be made much more realistic. The approximate theory of figures approach presented in this paper is 
readily generalized to arbitrary radial density profiles in the envelope, and the exact solution can also be 
extended to this case. 

6.1 Earth 

Table [3] presents the eccentricities and gravitational coefficient J2 for a two-layer model of the Earth with 
parameters pij p\ = 0.401, Qv = 0.1674, and 62 = 0.002. Results are given for the exact solution to the 
two-layer problem and for the theory of figures. Approximate solutions are valid to orders 1, 2, and 3 in 
the small parameter m. Table |3] also lists the observed values of _Ei, E2^ and J2. Two-layer models provide 
a good match to the observed eccentricities and an acceptable match to the gravitational coefficient. No 
attempt was made to fine tune the model parameters. For this case, even the theory of figures to first order 
in ra gives good agreement with the exact solution and with the observations. 

6.2 Mars 

Table [3] gives resuhs for a Mars model with pil px == 0.486, Qy = 0.125, £2 = 0.00347. There are no 
observations of the eccentricity of the Martian core-mantle boundary. The models provide good estimates 
of the eccentricity of the surface and the theory of figures approximations match the exact solution of E2 
quite closely. The eccentricity of the core-mantle boundary is less than that of the surface, as was the case 
for the Earth models, and E\ from the theory of figures approximations agrees rather well with the value of 
E\ from the exact solution. The model J2 is not in particularly good agreement with the observed Ji for 
Mars, but it is emphasized that we made no attempt to fine tune the model parameters to fit J2. Moreover, 
the radius of the Martian core and the densities of the Martian core and mantle are not known. 

6.3 Neptune 

Table[5]Hsts results for a Neptune model with pa/pi = 0.157334, Qy = 0.091125. £2 = 0.0254179 (parameter 
values based on a model in C.Z. Zhang (1997[)). It is emphasized that the two-layer model with a constant- 
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density envelope is not a good model for an ice giant planet like Neptune. Nevertheless, the shape of the 
surface is not too different from Neptune's observed shape, but J2 for the model is almost a factor of 2 larger 
than the observed value. We do not know if Neptune has a core-envelope configuration or a continuous radial 
density profile. 

6.4 Uranus 

Table l6] provides results for two Uranus models with model parameters given in the table and based on 



Horedt and Hubbard (19831. As was the case for Neptune, the Uranus models do okay in matching £'2 but 
fail to give good estimates for J2. Similar to Neptune, it is not known if the radial density profile of Uranus 
is a smooth one or if it contains a discontinuity associated with a core. 

7 Discussion and Conclusions 

The exact solution for the rotational distortion of a two-layer Maclaurin ellipsoid reported in [Kong et al.| 



(2010) has been extended here to provide formulas for the standard spherical harmonic expansion of the 
external gravitational field of the body. We have also presented a new approach to the evaluation of the 
theory of figures based on numerical integration of ordinary differential equations. 

The classical Radau-Darwin formula is a low order result from the theory of figures and its realm of 
validity has been evaluated for the two-layer model using the exact solution. It was found that the Radau- 
Darwin approximation is not valid for the rotational parameter 62 — r2^/(27rGp2) > 0.17 since the formula 
predicts a surface eccentricity that is smaller than the eccentricity of the core-envelope boundary. Interface 
eccentricity must be smaller than surface eccentricity. For an envelope density of 3000 kg m~'^ the failure of 
the Radau-Darwin formula corresponds to a rotation period of about 3 hr. 

The generalized Roche model, a two-layer model with an envelope density equal to zero, provides a simple 
model against which to evaluate the validity of the theory of figures against the exact solution. It was found 
that the theory of figures only slightly underestimates the eccentricities of the surface and core-envelope 
interface compared with the exact solution. 

Application of the exact solution and the theory of figures is made to models of Earth, Mars, Uranus, 
and Neptune. It is found that the two-layer model with constant densities in the layers can provide realistic 
approximations to terrestrial planets and icy outer planet satellites. This is perhaps not surprising since the 
zeroth order structure of these planetary bodies is similar to the two-layer model with constant densities in 
the layers. The situation is not as straightforward for giant planets since a constant density envelope is not 
a particularly good representation of the density in the outer layers of such planets. However, the theory of 
figures, as developed in this paper, is readily generalized to models with arbitrary radial density profiles in 
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the envelope (though we have not carried this out in this paper). Such models will be particularly useful for 
Jupiter and Saturn which might possess heavy element cores surrounded by gaseous envelopes. The envelope 
density can be represented by polynomial functions of radius. Inversions of gravitational data based on these 
models provide constraints on the gas giant interiors independent of assumptions about composition and 
equations of state. The exact solution for the two-layer Maclaurin ellipsoid can also be extended to allow 
for a non-constant radial profile of envelope density. This is not as straightforward as the generalization of 
the theory of figures, but it can be done. The solutions for two-layer bodies can therefore provide acceptable 
models for the rotational distortion of terrestrial, gas giant, and ice giant planetary bodies. These solutions 
can also serve as benchmarks to test the validity of complicated numerical models that invert gravitational 
and shape data to infer the interior structure of planets. 
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Figure 1: Surfaces of constant total (gravity) potential for cases (a) Qv = 0.5, p\l pi = 2, £2 = 0.18 and (b) 
Qv = 0.25, pxIpi = 2, £2 = 0.05. 
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Figure 2: The eccentricity of the total potential isosurfaces in Fig. [I] plotted as a function of radius. 
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First Order Function Fi 
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Figure 3: Function Fi{j3) for normalized envelope density 5e of 0.5 and core radius fie of 0.5. 
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Second Order Function K2 
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Figure 4: Function i^2(/3) for normalized envelope density Se of 0.5 and core radius /3c of 0.5. 



34 



Second Order Function F2 
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Figure 5: Function F2{fi) for normalized envelope density 5e of 0.5 and core radius Pc of 0.5. 
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Third Order Function H3 
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Figure 6: Function H^{j3) for normalized envelope density 5e of 0.5 and core radius j3c of 0.5. 
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Third Order Function K3 
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Figure 7: Function K3{/3) for normalized envelope density Se of 0.5 and core radius /3c of 0.5. 
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Third Order Function F3 
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Figure 8: Function F^^/S) for normalized envelope density 6e of 0.5 and core radius f3c of 0.5. 



38 



Table 1: Eccentricities of the interface Ei and the surface E2 as a function of the rotation parameter £2 for 



Qv — 0-5 and pi/ P2 — 2. The surface eccentricity based on the Radau-Darwin approximation is E,^ 



R-D 



£2 



El 



Eo 



E. 



R-D 



0.01000 


0.13959 


0.14390 


0.14383 


0.02000 


0.19761 


0.20330 


0.20288 


0.03000 


0.24158 


0.24860 


0.24782 


0.04000 


0.27886 


0.28670 


0.28540 


0.05000 


0.31147 


0.32010 


0.31824 


0.06000 


0.34050 


0.35010 


0.34769 


0.07000 


0.36770 


0.37770 


0.37453 


0.08000 


0.39250 


0.40320 


0.39931 


0.09000 


0.41592 


0.42710 


0.42238 


0.10000 


0.43790 


0.44960 


0.44402 


0.11000 


0.45864 


0.47090 


0.46441 


0.12000 


0.47895 


0.49130 


0.48372 


0.13000 


0.49789 


0.51070 


0.50208 


0.14000 


0.51604 


0.52930 


0.51957 


0.15000 


0.53399 


0.54730 


0.53630 


0.16000 


0.55103 


0.56460 


0.55232 


0.17000 


0.56739 


0.58130 


0.56771 


0.18000 


0.58335 


0.59750 


0.58250 


0.19000 


0.59873 


0.61320 


0.59674 


0.20000 


0.61401 


0.62854 


0.61047 


0.21000 


0.62897 


0.64350 


0.62373 


0.22000 


0.64315 


0.65800 


0.63654 


0.23000 


0.65732 


0.67224 


0.64894 


0.24000 


0.67128 


0.68620 


0.66093 
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Table 2: Comparison of interface and surface eccentricities for several models with small envelope densities 
computed from the exact theory and the theory of figures Roche model evaluated to second order in m. 



P2/P1 


<^2P2/Pl 


Qv 


El 


E2 


Ej°^P2=(^) 


i?J°F(p2-0) 


10-2 


0.05 


0.5 


0.4297 


0.4629 


0.4268 


0.4603 


10-3 


0.05 


0.5 


0.4292 


0.46296 


0.4268 


0.4603 


10-4 


0.05 


0.5 


0.4291 


0.46296 


0.4268 


0.4603 


10-2 


0.05 


0.33 


0.4378 


0.5189 


0.4268 


0.5141 


10-3 


0.05 


0.33 


0.4366 


0.5195 


0.4268 


0.5141 


10-4 


0.05 


0.33 


0.4365 


0.51955 


0.4268 


0.5141 


10-2 


0.02 


0.5 


0.2730 


0.2953 


0.2724 


0.2949 


10-3 


0.02 


0.5 


0.27297 


0.29537 


0.2724 


0.2949 


10-4 


0.02 


0.5 


0.2728 


0.29537 


0.2724 


0.2949 


10-2 


0.02 


0.33 


0.2757 


0.3314 


0.2724 


0.3313 


10-3 


0.02 


0.33 


0.2749 


0.3318 


0.2724 


0.3313 


10-4 


0.02 


0.33 


0.2748 


0.3318 


0.2724 


0.3313 
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Table 3: Application of exact solution and theory of figures to Earth. 
Two-layer Earth model: pa/pi == 0.401, Qy = 0.1674, £3 = 0.002. 





1st order 


2nd order 


3rd order 


exact 


observed 


El 


0.070765 


0.0707906 


0.0707907 


0.070593 


0.0707 


E2 


0.0810949 


0.0811186 


0.0811188 


0.081000 


0.082 


J2(106) 


1115.25 


1114.19 


1114.19 


1110.2 


1080 
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Table 4: Application of exact solution and theory of figures to Mars. 
Two-layer Mars model: pa/pi = 0.486, Qv = 0.125, ea = 0.00347. 





1st order 


2nd order 


3rd order 


exact 


observed 


El 


0.0888247 


0.0888743 


0.0888747 


0.088859 




E2 


0.100246 


0.100294 


0.100295 


0.10030 


0.10837 


J2(106) 


1825.82 


1823.18 


1823.18 


1823.1 


1960.0 
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Table 5: Application of exact solution and theory of figures to Neptune. Observed values of eccentricity and 
gravitational coefficient are from Lindal ( |1992 ) and Jacobson ( 2009 ) . 

Two-layer Neptune model: p^j px = 0.157334, Qy = 0.091125, £2 = 0.0254179. 





1st order 


2nd order 


3rd order 


exact 


observed 


El 


0.143134 


0.143506 


0.143515 


0.15147 


- 


E2 


0.209326 


0.209642 


0.209658 


0.21019 


0.18414 


J2{10') 


6228.69 


6188.61 


6188.92 


6241.0 


3408 
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Table 6: Application of exact solution and theory of figures to Uranus. Observed values of eccentricity and 



gravitational coefficient are from Lindal (1992) and Jacobson (2007). 



Two-layer Uranus model: pa/pi = 0.0883529, Qv = 0.421875, €2 = 0.103112. 



1st order 



2nd order 



3rd order 



exact 



observed 



El 


0.186917 


0.187284 


0.187296 


0.18752 


- 


E2 


0.207279 


0.207599 


0.207616 


0.20780 


0.212918 


J2(10«) 


4847.54 


4812.63 


4812.62 


4821.4 


3341 



Two-layer Uranus model: pz/pi == 0.0791231, Qv = 0.0563272, 63 = 0.0318902. 





1st order 


2nd order 


3rd order 


exact 


observed 


El 


0.115322 


0.115648 


0.115655 


0.14160 


- 


E2 


0.213329 


0.213629 


0.213648 


0.21473 


0.212918 


J2{10') 


5718.07 


5679.99 


5680.32 


5801.4 


3341 
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A Level-Surface Coefficients for the Spheroidal Functions /, k and 
h to Order Three 

The radial coordinate r, normalized to the mean radius s, can be written as follows as a truncated power 
series to order three in m. 



-1 + mf(^- A + m\ (^ - V + 4/") + ^ {mff (4 - 33/i' + 27^.^) 
+m^h (^ - 4m^ + 9m^ - 5m«) + {mff (11 - g^^ + 5^^ - ^m^) 



(101) 



This expansion for r is equivalent to ( 32 1 , but with the powers of m emphasized and stated explicitly. With 
k and h set equal to zero, it is the expansion for an ellipse with flattening mf, and with origin of coordinates 
at the center of the ellipse. 

With this function for r/s, it is straightforward to derive the coefficients CL for the potential functions 
Vi to arbitrary order by means of ( 45 ) , and by means of the procedure used to derive the coefficients Cj, in 



(42 1. For order one {i — 1), the result for the spheroidal functions is obtained to order 2 in the form 



cl~ 


"5^ + 35^ +35^ 


ch 


7-' 7 ■' 35 


ch 


36 402 48 
35'' 385'' 385 


ch 


77-' 77 



(102) 
and for order two, the coefficients to first order are 

2 21-' 
^4_i , 200 
^4-^+231^ 



50 
33' 



Gt = ^^ (103) 



For order three, there is only one non-zero coefficient to order one, C% which is equal to one 



The coefficients for the mass contribution exterior to the level surface at /3 follow from (46). There is 



only one non-zero coefficient for order zero, the coefficient Cq' which is equal to one. For order one the 
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coefRcients can be written as 



° 15 315 105 

2 2r 105 

^ 35-' 55^ 385 
Cf^'p' + ^^k (104) 



and for order 2 as, 



ct=o 



c^'=-^f 



C4 -1 " ^/ 



16 

21- 
160 

231- 
40 

33" 



C^'^-t;^! (105) 



and for order 3 there is just one non-zero coefficient Cg equal to one 



B Level-Surface Potential Functions 

The coefRcients C^^ and Cj'- derived in Appendix [aI can be substituted into (51 1 for the potential functions 
The internal normalized potential Aq on a level surface is obtained immediately to third order as 



8 .2 , 584 ,3 , 64 „,\ , /2 , , 13 
45-' 2835'' 315-' / V5 35^ 35 ' 



-S'„-l^f+^f' + ^k)S',+ i-, + ^f+^j' + ^k]n, (106) 



Similarly, the second degree potential function A2 is obtained immediately by (51), but it is simplified 
somewhat by multiplying it through by 3/2. Because it must be independent of /i on a level surface, and 
because it is multiplied by P2 (/i), it is equal to zero. The final expression for A2 is 

A^Jf+'^f + '-'f-h + h+^fAso 
V 42 63'' 7 7 105'' J 

Because A2 is zero, a solution for the small rotational parameter m can be found to third order. The 
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result is 



m=(2/-^/2 



49-' 7 



8 584 
7 ^ 245 



^/fc) So 



7'' 147-^ 7 



20 



/^4 



3 - 4/ + 2/2 



16, 



k]S'n 
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(108) 



For the higher degree potentials, the coefficients Cj, and C2J are substituted into (51 1. Then the above 



2i 



expression for m is substituted into the result, and terms higher than order three are dropped. In this way 
the centrifugal potential enters explicitly only in Aq and A2. When this procedure is applied to ^4 and the 
result is multiplied through by 35/4, the final expression is, 

2152 



A^ 
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— r h 

77 11 
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2385 2 156 
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'35 
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(109) 



Similarly for Aq, with the result multiplied through by —33/8, the final result is, 

32 
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25 



fSi 



33 



fk 5, 
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15,2^60 
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(110) 



Except for two obvious typographical errors in the m term for Aq, these expressions for Aq, A2, A4 and 



Aq agree with expressions given by Zharkov and Trubitsyn (1978). They can of course be carried to higher 



order, either by introducing higher-order spheroidal functions into ( 32 ) or by extending ( 38 1 to arbitrary 



order, as carried out by Zharkov and Trubitsyn to fifth order ( Zharkov and Trubitsyn 1978 ) 



C Evaluation of the Gravitational Moments 



The evaluation of the gravitational moments 821 and Sjj that appear in the potentials A2i is straightforward. 



An expression for r/s to arbitrary order is simply substituted into ( 53 1 and ( 54 ) and the integration is carried 



out over /i. The appropriate third-order expression for r/s is given by( 101 1, and the third-order expressions 



for the functions (pi and (j)'^ evaluate to the following. They agree with expressions given by Zharkov and 
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Trubitsyn ( 1978 ) 



1 = 1 
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(111) 
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D Functions Gji for the Differential Equations 



As described in section 5.2 the procedure for generating the ODE of (63) sequentially produces the right 



hand side of the equations as functions Gji of 13. The results of this process are 
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(112) 



These general equations, when applied to a particular problem, are not as complicated as they appear. The 



application of Eq. 112 to the two-zone model in Sec. 5.6 illustrates the method in more detail. It illustrates 
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our preferred method for application of the ODE approach to any interior calculation in general. 
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